library(tidyverse)
library(haven)
library(pals)
library(sandwich)
library(lmtest)
library(cowplot)
library(xtable)
library(flextable)

source("0_conjoint_functions.R")

#### -1. Plot formatting defaults ####

countries <- c("AUS","BR","CAN","CHL","CHN","COL","FR","GHA","IND","IT","JPN","SP","UGA","UK","US", "ZAF")
en_surveys <- c("AUS","CAN","IND","UGA","UK","US")
country_colours <- c("#ff0029","#377eb8",
                     "#66a61e", "#984ea3",
                     "#00d2d5", "#ff7f00",
                     "#af8d00", "#7f80cd",
                     "#b3e900", "#c42e60", 
                     "#a65628", "#f781bf", 
                     "#8dd3c7", "#bebada",
                     "#fb8072","#80b1d3")

plot_colours <- setNames(country_colours,
                         c("GHA","UGA","ZAF",
                           "BR","CAN","CHL","COL","US",
                           "AUS","CHN","IND","JPN",
                           "FR","IT","SP","UK"))

plot_theme <- theme_bw() + theme(
  axis.text = element_text(size = 11, color = "grey40"),
  axis.title = element_text(size = 11),
  legend.text = element_text(size = 11, color = "grey40")
)

#### 0. READ & FORMAT CONJOINT DATA ####
ref_levels = c(gdp = "0%",
               jobs = "0%",
               supplies = "Quick to procure",
               deaths = "10 per million",
               vaccinated = "5%",
               lockdown = "10 weeks")

# See format_data.R for data conversion/translation 

conjoint_data <- read_csv("formatted_data/conjoint_data.csv") %>% 
  set_refs() # set attribute reference levels

mod_ctrls <- c("gender","age","EDUCATION_LEVEL", "income_abovemed","marital_status","dep_children")

#### CONJOINT ANALYSIS ####

##### FIGURE 1 #####


## GLOBAL MODEL 

pooled_fe_mod <- conjoint_lm(data = conjoint_data, country_fe = TRUE, controls = mod_ctrls)

## RESULTS BY COUNTRY ##

country_mods <- lapply(countries, function (x) {
  
  country_mod <- conjoint_lm(conjoint_data[conjoint_data$country == x, ], 
                             country_fe = FALSE,
                             controls = mod_ctrls) %>% 
    format_model(., ref_levels, cluster = "c_id", label = x) %>% 
    mutate(att = case_when(att == "GDP growth" ~ "GDP",
                           att == "Job growth" ~ "Jobs",
                           att == "Lockdown length" ~ "Lockdown",
                           att == "Vaccination rate" ~ "Vaccinated",
                           att == "Vaccine procurement" ~ "Supplies",
                           !is.na(att) ~ "Deaths"))
}) %>% do.call("rbind", .) %>% 
  
  # Add in global results
  rbind(format_model(pooled_fe_mod, ref_levels, cluster = "country", label= "Pooled") %>% 
          mutate(att = case_when(att == "GDP growth" ~ "GDP",
                                 att == "Job growth" ~ "Jobs",
                                 att == "Lockdown length" ~ "Lockdown",
                                 att == "Vaccination rate" ~ "Vaccinated",
                                 att == "Vaccine procurement" ~ "Supplies",
                                 !is.na(att) ~ "Deaths"))) %>% 
  # format results
  mutate(region = case_when(model %in% c("GHA", "UGA", "ZAF") ~ "Africa",
                            model %in% c("FR","IT","UK","SP") ~ "Europe",
                            model %in% c("AUS","CHN","IND", "JPN") ~ "APac",
                            model %in% c("BR","CAN","CHL","COL","US") ~ "Americas",
                            model == "Pooled" ~ "Global"),
         model_lab = ifelse(non_ref==1, model, "Reference")) %>% 
  mutate(model_lab = relevel(as.factor(model_lab), ref = "Reference"))

region_data <- split(country_mods,
                     country_mods$region)

fig_limits <- country_mods %>% 
  mutate(x_low = est - 1.96*std_error,
         x_high = est + 1.96*std_error) %>% 
  summarise(x_low = min(x_low, na.rm = TRUE)-0.01,
            x_high = max(x_high, na.rm = TRUE) + 0.01)

regions <- str_sort(unique(country_mods$region))

conjoint_colours <- c(plot_colours, c("Pooled" = "#30c7d4"))

region_data <- region_data[c("Global","Africa","Americas","APac","Europe")]

region_plots <- lapply(1:length(region_data), function (i) {
  
  region_data[[i]] %>% 
    arrange(model_lab) %>% 
    ggplot(aes(x = est, y = coef, xmin = lower95, xmax = upper95)) +
    geom_point(position = position_dodge(width = -1), size = 1.5, alpha = 0.7) +
    geom_errorbarh(height = 0, position = position_dodge(width = -1), size = 0.8, alpha = 0.7) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    facet_grid(att~region, space = "free", scales = "free_y") +
    aes(color = model_lab) +
    scale_color_manual(values=conjoint_colours) +
    labs(x = "AMCE", y = "", color = "") +
    scale_y_discrete(limits = rev(levels("term"))) +
    guides(color=guide_legend(ncol = 1, byrow = TRUE)) +
    xlim(fig_limits$x_low, fig_limits$x_high) +
    theme_bw()+
    theme(legend.position = "bottom",
          text = element_text(size = 11),
          legend.text = element_text(size = 11),
          strip.text = element_text(size = 11),
          plot.margin = unit(c(0, 0, 0, 0), "cm")) +
    
    if (i == 1) {
      
      theme(axis.text.y = element_text(size = 11),
            strip.background.y = element_blank(),
            strip.text.y = element_blank())
      
    } else if (i != length(regions)) {
      
      theme(axis.text.y = element_blank(),
            axis.ticks.y = element_blank(),
            strip.background.y = element_blank(),
            strip.text.y = element_blank(),
            plot.margin=margin(l=-0.4,unit="cm"))
      
    } else {
      
      theme(axis.text.y = element_blank(),
            axis.ticks.y = element_blank(),
            plot.margin=margin(l=-0.4,unit="cm"))
      
    }
}
)

country_plot <- plot_grid(plotlist = region_plots, nrow = 1,
          align = "h", rel_widths = c(1.9,1,1,1,1.1))

plot(country_plot)

ggsave(country_plot, filename = "figures/figure1.pdf", width = 9, height = 11)

##### SUBGROUP ANALYSIS #####
###### FIGURE S1 -- CONJOINT SCREENSHOT #####
## Produced by hand
###### FIGURE S2 -- RESULTS BY AGE ###### 

age_models <- conjoint_data %>% 
  mutate(age_bin = ifelse(age >= 50, "50+ years old", "Less than 50 years old")) %>% 
  group_by(age_bin) %>% 
  nest() %>% 
  filter(!is.na(age_bin)) %>% 
  mutate(model = map(data, ~conjoint_lm(data = .x, country_fe = TRUE, controls = mod_ctrls[mod_ctrls!="age"])),
         model_table = map2(model, age_bin, ~format_model(.x, ref_levels, label = .y, cluster = "country")))

age_results <- do.call("rbind", age_models$model_table)

age_plot <- ggplot(age_results, aes(x = est, y = coef, xmin = lower95, xmax = upper95, color = model, alpha = non_ref)) +
  geom_point(size = 2.5, position = position_dodge(width=0.8)) +
  geom_errorbarh(height = 0, size = 0.8, position = position_dodge(width=0.8)) +
  facet_wrap(att ~ ., ncol = 1, scales = "free_y", strip.position = "top") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_alpha_manual(values = c(0.3,1), guide = "none") +
  theme_bw() +
  theme(legend.position = "bottom") + 
  labs(x = "AMCE",y = "")

plot(age_plot)
ggsave(age_plot, filename = "figures/age_plot.pdf", width = 8, height = 7)


###### FIGURE S3 -- RESULTS BY EDUCATION LEVEL ###### 
education_mods <- conjoint_data %>% 
  filter(EDUCATION_LEVEL != "Unknown/missing") %>% 
  mutate(education = ifelse(EDUCATION_LEVEL %in% c("Less than primary completed","Primary completed"),
                            "Up to primary completed",
                            EDUCATION_LEVEL),
         education = relevel(factor(education), ref = "Up to primary completed")) %>% 
  group_by(education) %>% 
  nest() %>% 
  mutate(model = map(data, ~conjoint_lm(data = .x, country_fe = TRUE, controls = mod_ctrls[mod_ctrls!="EDUCATION_LEVEL"])),
         model_table = map2(model, education, ~format_model(.x, ref_levels, label = .y, cluster = "country")))

education_results <- do.call("rbind", education_mods$model_table)

education_plot <- ggplot(education_results, 
                         aes(x = est, y = coef, xmin = lower95, xmax = upper95, 
                             color = model, alpha = non_ref)) +
  geom_point(size = 2.5, position = position_dodge(width=0.8)) +
  geom_errorbarh(height = 0, size = 0.8, position = position_dodge(width=0.8)) +
  facet_wrap(att ~ ., ncol = 1, scales = "free_y", strip.position = "top") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_alpha_manual(values = c(0.3,1), guide = "none") +
  theme_bw() +
  theme(legend.position = "bottom") + 
  labs(x = "AMCE",y = "", color = "Education level")

plot(education_plot)
ggsave(education_plot, filename = "figures/education_level_plot.pdf", width = 8, height = 7)

###### FIGURE S4 -- RESULTS BY GENDER ###### 
gender_mods <- conjoint_data %>% 
  filter(gender %in% c("Male","Female")) %>% 
  group_by(gender) %>% 
  nest() %>% 
  mutate(model = map(data, ~conjoint_lm(data = .x, country_fe = TRUE, controls = mod_ctrls[mod_ctrls!="gender"])),
         model_table = map2(model, gender, ~format_model(.x, ref_levels, label = .y, cluster = "country")))

gender_results <- do.call("rbind", gender_mods$model_table)

gender_plot <- ggplot(gender_results, 
                         aes(x = est, y = coef, xmin = lower95, xmax = upper95, 
                             color = model, alpha = non_ref)) +
  geom_point(size = 2.5, position = position_dodge(width=0.8)) +
  geom_errorbarh(height = 0, size = 0.8, position = position_dodge(width=0.8)) +
  facet_wrap(att ~ ., ncol = 1, scales = "free_y", strip.position = "top") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_alpha_manual(values = c(0.3,1), guide = "none") +
  theme_bw() +
  theme(legend.position = "bottom") + 
  labs(x = "AMCE",y = "", color = "Education level")

plot(gender_plot)
ggsave(gender_plot, filename = "figures/gender_plot.pdf", width = 8, height = 7)

###### FIGURE S5 -- RESULTS BY COVID EXPOSURE (RETROSPECTIVE) ###### 

exp_models <- conjoint_data %>% 
  mutate(covid_exp = case_when(vac_hes_4_1 == "Yes" ~ "Infected in last year", 
                               vac_hes_4_1 == "No" ~ "No infection")) %>% 
  group_by(covid_exp) %>% 
  nest() %>% 
  filter(!is.na(covid_exp)) %>% 
  mutate(model = map(data, ~conjoint_lm(data = .x, country_fe = TRUE, controls = mod_ctrls)),
         model_table = map2(model, covid_exp, ~format_model(.x, ref_levels, label = .y, cluster = "country")))

exp_results <- do.call("rbind", exp_models$model_table)

exp_plot <- ggplot(exp_results, 
                   aes(x = est, y = coef, xmin = lower95, xmax = upper95, 
                       color = model, alpha = non_ref)) +
  geom_point(size = 2.5, position = position_dodge(width=0.8)) +
  geom_errorbarh(height = 0, size = 0.8, position = position_dodge(width=0.8)) +
  facet_wrap(att ~ ., ncol = 1, scales = "free_y", strip.position = "top") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_alpha_manual(values = c(0.3,1), guide = "none") +
  theme_bw() +
  theme(legend.position = "bottom") + 
  labs(x = "AMCE",y = "", color = "Previous Covid-19 infection")

plot(exp_plot)
ggsave(exp_plot, filename = "figures/covid_exposure_plot.pdf", width = 8, height = 7)

###### FIGURE S6 -- RESULTS BY COVID EXPOSURE (PROSPECTIVE) ###### 

exp_prospective <- conjoint_data %>% 
  mutate(covid_exp_p = case_when(vac_hes_7 > 50 ~ "Expects to be infected", 
                               vac_hes_7 <= 50 ~ "No expectation of infection")) %>% 
  group_by(covid_exp_p) %>% 
  nest() %>% 
  filter(!is.na(covid_exp_p)) %>% 
  mutate(model = map(data, ~conjoint_lm(data = .x, country_fe = TRUE, controls = mod_ctrls)),
         model_table = map2(model, covid_exp_p, ~format_model(.x, ref_levels, label = .y, cluster = "country")))

exp_p_results <- do.call("rbind", exp_prospective$model_table)

exp_p_plot <- ggplot(exp_p_results, 
                   aes(x = est, y = coef, xmin = lower95, xmax = upper95, 
                       color = model, alpha = non_ref)) +
  geom_point(size = 2.5, position = position_dodge(width=0.8)) +
  geom_errorbarh(height = 0, size = 0.8, position = position_dodge(width=0.8)) +
  facet_wrap(att ~ ., ncol = 1, scales = "free_y", strip.position = "top") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_alpha_manual(values = c(0.3,1), guide = "none") +
  theme_bw() +
  theme(legend.position = "bottom") + 
  labs(x = "AMCE",y = "", color = "Prospective judgement")

plot(exp_p_plot)
ggsave(exp_p_plot, filename = "figures/covid_prospective_exposure_plot.pdf", width = 8, height = 7)

###### FIGURE S7 -- RESULTS BY COVID CASES ###### 

covid_cases <- conjoint_data %>% 
  mutate(cases_w2_p = cases_w2/pop2019_w2,
         cases_bin = cases_w2_p > median(cases_w2_p, na.rm = TRUE)) %>% 
  group_by(cases_bin) %>% 
  nest() %>% 
  filter(!is.na(cases_bin)) %>% 
  mutate(model = map(data, ~conjoint_lm(data = .x, country_fe = TRUE, controls = mod_ctrls)),
         model_table = map2(model, cases_bin, ~format_model(.x, ref_levels, label = .y, cluster = "country")))

covid_cases_results <- do.call("rbind", covid_cases$model_table)

covid_cases_plot <- ggplot(covid_cases_results, 
                     aes(x = est, y = coef, xmin = lower95, xmax = upper95, 
                         color = model, alpha = non_ref)) +
  geom_point(size = 2.5, position = position_dodge(width=0.8)) +
  geom_errorbarh(height = 0, size = 0.8, position = position_dodge(width=0.8)) +
  facet_wrap(att ~ ., ncol = 1, scales = "free_y", strip.position = "top") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_alpha_manual(values = c(0.3,1), guide = "none") +
  theme_bw() +
  theme(legend.position = "bottom") + 
  labs(x = "AMCE",y = "", color = "Cases per capita above median?")

plot(covid_cases_plot)
ggsave(covid_cases_plot, filename = "figures/covid_cases_plot.pdf", width = 8, height = 7)


###### FIGURE S8 -- RESULTS BY COVID DEATHS ###### 

covid_deaths <- conjoint_data %>% 
  mutate(deaths_w2_p = deaths_w2/pop2019_w2,
         deaths_bin = deaths_w2_p > median(deaths_w2_p, na.rm = TRUE)) %>% 
  group_by(deaths_bin) %>% 
  nest() %>% 
  filter(!is.na(deaths_bin)) %>% 
  mutate(model = map(data, ~conjoint_lm(data = .x, country_fe = TRUE, controls = mod_ctrls)),
         model_table = map2(model, deaths_bin, ~format_model(.x, ref_levels, label = .y, cluster = "country")))

covid_deaths_results <- do.call("rbind", covid_deaths$model_table)

covid_deaths_plot <- ggplot(covid_deaths_results, 
                           aes(x = est, y = coef, xmin = lower95, xmax = upper95, 
                               color = model, alpha = non_ref)) +
  geom_point(size = 2.5, position = position_dodge(width=0.8)) +
  geom_errorbarh(height = 0, size = 0.8, position = position_dodge(width=0.8)) +
  facet_wrap(att ~ ., ncol = 1, scales = "free_y", strip.position = "top") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_alpha_manual(values = c(0.3,1), guide = "none") +
  theme_bw() +
  theme(legend.position = "bottom") + 
  labs(x = "AMCE",y = "", color = "Deaths per capita above median?")

plot(covid_deaths_plot)
ggsave(covid_deaths_plot, filename = "figures/covid_deaths_plot.pdf", width = 8, height = 7)

###### FIGURE S9 -- RESULTS BY FOOD INDEX ###### 

food_mean <- mean(conjoint_data$food_indx, na.rm = TRUE)

food <- conjoint_data %>% 
  mutate(food_bin = case_when(food_indx > food_mean ~ "High impact on food resources", 
                              food_indx <= food_mean ~ "Low impact on food resources")) %>% 
  group_by(food_bin) %>% 
  nest() %>% 
  filter(!is.na(food_bin)) %>% 
  mutate(model = map(data, ~conjoint_lm(data = .x, country_fe = TRUE, controls = mod_ctrls)),
         model_table = map2(model, food_bin, ~format_model(.x, ref_levels, label = .y, cluster = "country")))

food_results <- do.call("rbind", food$model_table)

food_plot <- ggplot(food_results, 
                     aes(x = est, y = coef, xmin = lower95, xmax = upper95, 
                         color = model, alpha = non_ref)) +
  geom_point(size = 2.5, position = position_dodge(width=0.8)) +
  geom_errorbarh(height = 0, size = 0.8, position = position_dodge(width=0.8)) +
  facet_wrap(att ~ ., ncol = 1, scales = "free_y", strip.position = "top") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_alpha_manual(values = c(0.3,1), guide = "none") +
  theme_bw() +
  theme(legend.position = "bottom") + 
  labs(x = "AMCE",y = "", color = "Impact")

plot(food_plot)
ggsave(food_plot, filename = "figures/food_plot.pdf", width = 8, height = 7)

###### FIGURE S10 -- RESULTS BY FOOD PCA ###### 

food_pca <- conjoint_data %>% 
  mutate(food_pca = case_when(food_indx > 0 ~ "High impact on food resources", 
                              food_indx <= 0 ~ "Low impact on food resources")) %>% 
  group_by(food_pca) %>% 
  nest() %>% 
  filter(!is.na(food_pca)) %>% 
  mutate(model = map(data, ~conjoint_lm(data = .x, country_fe = TRUE, controls = mod_ctrls)),
         model_table = map2(model, food_pca, ~format_model(.x, ref_levels, label = .y, cluster = "country")))

food_pca_results <- do.call("rbind", food_pca$model_table)

food_pca_plot <- ggplot(food_pca_results, 
                    aes(x = est, y = coef, xmin = lower95, xmax = upper95, 
                        color = model, alpha = non_ref)) +
  geom_point(size = 2.5, position = position_dodge(width=0.8)) +
  geom_errorbarh(height = 0, size = 0.8, position = position_dodge(width=0.8)) +
  facet_wrap(att ~ ., ncol = 1, scales = "free_y", strip.position = "top") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_alpha_manual(values = c(0.3,1), guide = "none") +
  theme_bw() +
  theme(legend.position = "bottom") + 
  labs(x = "AMCE",y = "", color = "Impact")

plot(food_pca_plot)
ggsave(food_pca_plot, filename = "figures/food_pca_plot.pdf", width = 8, height = 7)

###### FIGURE S11 -- RESULTS BY IDEOLOGY ###### 

ideology <- conjoint_data %>% 
  mutate(ideo_bin = case_when(ideology > 5 ~ "Right", 
                              ideology <= 5 ~ "Left")) %>% 
  group_by(ideo_bin) %>% 
  nest() %>% 
  filter(!is.na(ideo_bin)) %>% 
  mutate(model = map(data, ~conjoint_lm(data = .x, country_fe = TRUE, controls = mod_ctrls)),
         model_table = map2(model, ideo_bin, ~format_model(.x, ref_levels, label = .y, cluster = "country")))

ideo_results <- do.call("rbind", ideology$model_table)

ideo_plot <- ggplot(ideo_results, 
                        aes(x = est, y = coef, xmin = lower95, xmax = upper95, 
                            color = model, alpha = non_ref)) +
  geom_point(size = 2.5, position = position_dodge(width=0.8)) +
  geom_errorbarh(height = 0, size = 0.8, position = position_dodge(width=0.8)) +
  facet_wrap(att ~ ., ncol = 1, scales = "free_y", strip.position = "top") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_alpha_manual(values = c(0.3,1), guide = "none") +
  theme_bw() +
  theme(legend.position = "bottom") + 
  labs(x = "AMCE",y = "", color = "Ideology")

plot(ideo_plot)
ggsave(ideo_plot, filename = "figures/ideo_plot.pdf", width = 8, height = 7)

###### FIGURE S12 -- RESULTS BY HH INC ABOVE MEDIAN ###### 

hh_inc_abovemed <- conjoint_data %>% 
  group_by(income_abovemed) %>% 
  nest() %>% 
  filter(!is.na(income_abovemed)) %>% 
  mutate(model = map(data, ~conjoint_lm(data = .x, country_fe = TRUE, controls = mod_ctrls[mod_ctrls!="income_abovemed"])),
         model_table = map2(model, income_abovemed, ~format_model(.x, ref_levels, label = .y, cluster = "country")))

hh_inc_abovemed_results <- do.call("rbind", hh_inc_abovemed$model_table) %>% 
  mutate(model = ifelse(model == 1, "Above country median","Below country median"))
  
hh_inc_abovemed_plot <- ggplot(hh_inc_abovemed_results, 
                            aes(x = est, y = coef, xmin = lower95, xmax = upper95, 
                                color = as.factor(model), alpha = non_ref)) +
  geom_point(size = 2.5, position = position_dodge(width=0.8)) +
  geom_errorbarh(height = 0, size = 0.8, position = position_dodge(width=0.8)) +
  facet_wrap(att ~ ., ncol = 1, scales = "free_y", strip.position = "top") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_alpha_manual(values = c(0.3,1), guide = "none") +
  theme_bw() +
  theme(legend.position = "bottom") + 
  labs(x = "AMCE",y = "", color = "HH Income")

plot(hh_inc_abovemed_plot)
ggsave(hh_inc_abovemed_plot, filename = "figures/hh_inc_abovemed_plot.pdf", width = 8, height = 7)

###### FIGURE S13 -- RESULTS BY HH INCOME CHANGE ###### 

hh_inc_delta <- conjoint_data %>% 
  mutate(hh_inc_d = ifelse(hh_inc_delta %in% c("Do not know","Prefer not to say",NA), NA, hh_inc_delta)) %>% 
  group_by(hh_inc_d) %>% 
  nest() %>% 
  filter(!is.na(hh_inc_d)) %>% 
  mutate(model = map(data, ~conjoint_lm(data = .x, country_fe = TRUE,controls = mod_ctrls)),
         model_table = map2(model, hh_inc_d, ~format_model(.x, ref_levels, label = .y, cluster = "country")))

hh_inc_delta_results <- do.call("rbind", hh_inc_delta$model_table)

hh_inc_delta_plot <- ggplot(hh_inc_delta_results, 
                    aes(x = est, y = coef, xmin = lower95, xmax = upper95, 
                        color = model, alpha = non_ref)) +
  geom_point(size = 2.5, position = position_dodge(width=0.8)) +
  geom_errorbarh(height = 0, size = 0.8, position = position_dodge(width=0.8)) +
  facet_wrap(att ~ ., ncol = 1, scales = "free_y", strip.position = "top") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_alpha_manual(values = c(0.3,1), guide = "none") +
  theme_bw() +
  theme(legend.position = "bottom") + 
  labs(x = "AMCE",y = "", color = "Change in household income")

plot(hh_inc_delta_plot)
ggsave(hh_inc_delta_plot, filename = "figures/hh_inc_delta_plot.pdf", width = 8, height = 7)


###### FIGURE S14 -- RESULTS BY EQ5D VAS (TODAY) ###### 

eq5d_today <- conjoint_data %>% 
  mutate(eq5d_rate_now_bin = ifelse(eq5d_rate_now > 50, "51-100", "0-50")) %>% 
  mutate(eq5d_rate_now_bin = factor(eq5d_rate_now_bin, levels = c("0-50","51-100"))) %>% 
  group_by(eq5d_rate_now_bin) %>% 
  nest() %>% 
  filter(!is.na(eq5d_rate_now_bin)) %>% 
  mutate(model = map(data, ~conjoint_lm(data = .x, country_fe = TRUE,controls = mod_ctrls)),
         model_table = map2(model, eq5d_rate_now_bin, ~format_model(.x, ref_levels, label = .y, cluster = "country")))

eq5d_today_results <- do.call("rbind", eq5d_today$model_table)

eq5d_today_plot <- ggplot(eq5d_today_results, 
                          aes(x = est, y = coef, xmin = lower95, xmax = upper95, 
                              color = model, alpha = non_ref)) +
  geom_point(size = 2.5, position = position_dodge(width=0.8)) +
  geom_errorbarh(height = 0, size = 0.8, position = position_dodge(width=0.8)) +
  facet_wrap(att ~ ., ncol = 1, scales = "free_y", strip.position = "top") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_alpha_manual(values = c(0.3,1), guide = "none") +
  theme_bw() +
  theme(legend.position = "bottom") + 
  labs(x = "AMCE",y = "", color = "EQ5D VAS (Today)")

plot(eq5d_today_plot)
ggsave(eq5d_today_plot, filename = "figures/eq5d_today_plot.pdf", width = 8, height = 7)


###### FIGURE S15 -- RESULTS BY EQ5D VAS (DELTA) ###### 

eq5d_delta <- conjoint_data %>% 
  mutate(eq5d_delta_bin = ifelse(eq5d_rate_delta > 0, "Health improved", "Health the same/worse")) %>% 
  group_by(eq5d_delta_bin) %>% 
  nest() %>% 
  filter(!is.na(eq5d_delta_bin)) %>% 
  mutate(model = map(data, ~conjoint_lm(data = .x, country_fe = TRUE,controls = mod_ctrls)),
         model_table = map2(model, eq5d_delta_bin, ~format_model(.x, ref_levels, label = .y, cluster = "country")))

eq5d_delta_results <- do.call("rbind", eq5d_delta$model_table)

eq5d_delta_plot <- ggplot(eq5d_delta_results, 
                          aes(x = est, y = coef, xmin = lower95, xmax = upper95, 
                              color = model, alpha = non_ref)) +
  geom_point(size = 2.5, position = position_dodge(width=0.8)) +
  geom_errorbarh(height = 0, size = 0.8, position = position_dodge(width=0.8)) +
  facet_wrap(att ~ ., ncol = 1, scales = "free_y", strip.position = "top") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_alpha_manual(values = c(0.3,1), guide = "none") +
  theme_bw() +
  theme(legend.position = "bottom") + 
  labs(x = "AMCE",y = "", color = "EQ5D VAS (Change over 1 year)")

plot(eq5d_delta_plot)
ggsave(eq5d_delta_plot, filename = "figures/eq5d_delta_plot.pdf", width = 8, height = 7)

###### FIGURE S16 -- RESULTS BY WHO INDEX ###### 

who_mean <- mean(conjoint_data$who_indx, na.rm = TRUE)

who_indx <- conjoint_data %>% 
  mutate(who_indx_bin = case_when(who_indx > who_mean ~ "Above average health", 
                                  who_indx <= who_mean ~ "Below average health")) %>% 
  group_by(who_indx_bin) %>% 
  nest() %>% 
  filter(!is.na(who_indx_bin)) %>% 
  mutate(model = map(data, ~conjoint_lm(data = .x, country_fe = TRUE, controls = mod_ctrls)),
         model_table = map2(model, who_indx_bin, ~format_model(.x, ref_levels, label = .y, cluster = "country")))

who_indx_results <- do.call("rbind", who_indx$model_table)

who_indx_plot <- ggplot(who_indx_results, 
                    aes(x = est, y = coef, xmin = lower95, xmax = upper95, 
                        color = model, alpha = non_ref)) +
  geom_point(size = 2.5, position = position_dodge(width=0.8)) +
  geom_errorbarh(height = 0, size = 0.8, position = position_dodge(width=0.8)) +
  facet_wrap(att ~ ., ncol = 1, scales = "free_y", strip.position = "top") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_alpha_manual(values = c(0.3,1), guide = "none") +
  theme_bw() +
  theme(legend.position = "bottom") + 
  labs(x = "AMCE",y = "", color = "Self-reported WHO health assessment")

plot(who_indx_plot)
ggsave(who_indx_plot, filename = "figures/who_indx_plot.pdf", width = 8, height = 7)

###### FIGURE S17 -- RESULTS BY WHO PCA ###### 

who_pca <- conjoint_data %>% 
  mutate(who_pca_bin = case_when(who_pca > 0 ~ "Above average health score (PCA)", 
                              who_pca <= 0 ~ "Below average health score (PCA)")) %>% 
  group_by(who_pca_bin) %>% 
  nest() %>% 
  filter(!is.na(who_pca_bin)) %>% 
  mutate(model = map(data, ~conjoint_lm(data = .x, country_fe = TRUE,controls = mod_ctrls)),
         model_table = map2(model, who_pca_bin, ~format_model(.x, ref_levels, label = .y, cluster = "country")))

who_pca_results <- do.call("rbind", who_pca$model_table)

who_pca_plot <- ggplot(who_pca_results, 
                        aes(x = est, y = coef, xmin = lower95, xmax = upper95, 
                            color = model, alpha = non_ref)) +
  geom_point(size = 2.5, position = position_dodge(width=0.8)) +
  geom_errorbarh(height = 0, size = 0.8, position = position_dodge(width=0.8)) +
  facet_wrap(att ~ ., ncol = 1, scales = "free_y", strip.position = "top") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_alpha_manual(values = c(0.3,1), guide = "none") +
  theme_bw() +
  theme(legend.position = "bottom") + 
  labs(x = "AMCE",y = "", color = "WHO score (PCA)")

plot(who_pca_plot)
ggsave(who_pca_plot, filename = "figures/who_pca_plot.pdf", width = 8, height = 7)

###### FIGURE S18 -- RESULTS BY VACCINATION STATUS ###### 

vacc_mods <- conjoint_data %>%
  group_by(subj_vaccinated) %>% 
  nest() %>% 
  filter(!is.na(subj_vaccinated)) %>% 
  mutate(model = map(data, ~conjoint_lm(data = .x, country_fe = TRUE,controls = mod_ctrls)),
         model_table = map2(model, subj_vaccinated, ~format_model(.x, ref_levels, label = .y, cluster = "country")))

vacc_results <- do.call("rbind", vacc_mods$model_table) %>% 
  mutate(model = ifelse(model == 1, "Yes","No"))

vacc_plot <- ggplot(vacc_results, 
                       aes(x = est, y = coef, xmin = lower95, xmax = upper95, 
                           color = model, alpha = non_ref)) +
  geom_point(size = 2.5, position = position_dodge(width=0.8)) +
  geom_errorbarh(height = 0, size = 0.8, position = position_dodge(width=0.8)) +
  facet_wrap(att ~ ., ncol = 1, scales = "free_y", strip.position = "top") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_alpha_manual(values = c(0.3,1), guide = "none") +
  theme_bw() +
  theme(legend.position = "bottom") + 
  labs(x = "AMCE",y = "", color = "Vaccinated?")

plot(vacc_plot)

ggsave(vacc_plot, filename = "figures/vacc_plot.pdf", width = 8, height = 7)

###### FIGURE S19 -- RESULTS BY REFUSED VACCINE ###### 

vacc_refuse_mods <- conjoint_data %>%
  group_by(subj_vacc_refused) %>% 
  nest() %>% 
  filter(!is.na(subj_vacc_refused)) %>% 
  mutate(model = map(data, ~conjoint_lm(data = .x, country_fe = TRUE,controls = mod_ctrls)),
         model_table = map2(model, subj_vacc_refused, ~format_model(.x, ref_levels, label = .y, cluster = "country")))

vacc_refuse_results <- do.call("rbind", vacc_refuse_mods$model_table) %>% 
  mutate(model = ifelse(model == 1, "Yes","No"))

vacc_refuse_plot <- ggplot(vacc_refuse_results, 
                    aes(x = est, y = coef, xmin = lower95, xmax = upper95, 
                        color = model, alpha = non_ref)) +
  geom_point(size = 2.5, position = position_dodge(width=0.8)) +
  geom_errorbarh(height = 0, size = 0.8, position = position_dodge(width=0.8)) +
  facet_wrap(att ~ ., ncol = 1, scales = "free_y", strip.position = "top") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_alpha_manual(values = c(0.3,1), guide = "none") +
  theme_bw() +
  theme(legend.position = "bottom") + 
  labs(x = "AMCE",y = "", color = "Refused vaccine?")

plot(vacc_refuse_plot)

ggsave(vacc_refuse_plot, filename = "figures/vacc_refuse_plot.pdf", width = 8, height = 7)

###### FIGURE S20 -- RESULTS BY CONTAINMENT POLICIES ###### 

containment_mods <- conjoint_data %>%
  mutate(containment_bin = ifelse(C_pca <= median(C_pca), "Low","High")) %>% 
  group_by(containment_bin) %>% 
  nest() %>% 
  filter(!is.na(containment_bin)) %>% 
  mutate(model = map(data, ~conjoint_lm(data = .x, country_fe = TRUE,controls = mod_ctrls)),
         model_table = map2(model, containment_bin, ~format_model(.x, ref_levels, label = .y, cluster = "country")))

containment_results <- do.call("rbind", containment_mods$model_table)

containment_plot <- ggplot(containment_results, 
                           aes(x = est, y = coef, xmin = lower95, xmax = upper95, 
                               color = model, alpha = non_ref)) +
  geom_point(size = 2.5, position = position_dodge(width=0.8)) +
  geom_errorbarh(height = 0, size = 0.8, position = position_dodge(width=0.8)) +
  facet_wrap(att ~ ., ncol = 1, scales = "free_y", strip.position = "top") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_alpha_manual(values = c(0.3,1), guide = "none") +
  theme_bw() +
  theme(legend.position = "bottom") + 
  labs(x = "AMCE",y = "", color = "Containment policies")

plot(containment_plot)

ggsave(containment_plot, filename = "figures/containment_plot.pdf", width = 8, height = 7)

###### FIGURE S21 -- RESULTS BY ECONOMIC POLICIES ###### 

economic_mods <- conjoint_data %>%
  mutate(economic_bin = ifelse(E_pca <= median(E_pca), "Low","High")) %>% 
  group_by(economic_bin) %>% 
  nest() %>% 
  filter(!is.na(economic_bin)) %>% 
  mutate(model = map(data, ~conjoint_lm(data = .x, country_fe = TRUE,controls = mod_ctrls)),
         model_table = map2(model, economic_bin, ~format_model(.x, ref_levels, label = .y, cluster = "country")))

economic_results <- do.call("rbind", economic_mods$model_table)

economic_plot <- ggplot(economic_results, 
                           aes(x = est, y = coef, xmin = lower95, xmax = upper95, 
                               color = model, alpha = non_ref)) +
  geom_point(size = 2.5, position = position_dodge(width=0.8)) +
  geom_errorbarh(height = 0, size = 0.8, position = position_dodge(width=0.8)) +
  facet_wrap(att ~ ., ncol = 1, scales = "free_y", strip.position = "top") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_alpha_manual(values = c(0.3,1), guide = "none") +
  theme_bw() +
  theme(legend.position = "bottom") + 
  labs(x = "AMCE",y = "", color = "Economic policies")

plot(economic_plot)

ggsave(economic_plot, filename = "figures/economic_plot.pdf", width = 8, height = 7)

###### FIGURE S31 -- RESULTS BY ATTENTION ######

attention_models <- conjoint_data %>% 
  mutate(attention = ifelse(attention == 1, "Yes", "No")) %>% 
  group_by(attention) %>% 
  nest() %>% 
  mutate(model = map(data, ~conjoint_lm(data = .x, country_fe = TRUE, controls = mod_ctrls)),
         model_table = map2(model, attention, ~format_model(.x, ref_levels, label = .y, cluster = "country")))

attention_results <- do.call("rbind", attention_models$model_table)

attention_plot <- ggplot(attention_results, aes(x = est, y = coef, xmin = lower95, xmax = upper95, color = model, alpha = non_ref)) +
  geom_point(size = 2.5, position = position_dodge(width=0.8)) +
  geom_errorbarh(height = 0, size = 0.8, position = position_dodge(width=0.8)) +
  facet_wrap(att ~ ., ncol = 1, scales = "free_y", strip.position = "top") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_alpha_manual(values = c(0.3,1), guide = "none") +
  theme_bw() +
  theme(legend.position = "bottom") + 
  labs(x = "AMCE",y = "", color = "Passed attention check?")

plot(attention_plot)
ggsave(attention_plot, filename = "figures/attention_plot.pdf", width = 8, height = 7)



##### FIGURE 2 #####

fig2_results <- rbind(
  food$model %>% 
    wald_test_submodel(.) %>% 
    filter(str_detect(coef, "vaccinated")) %>% 
    mutate(subgroup = "Food Poverty"),
  
  containment_mods$model %>% 
    wald_test_submodel(.) %>% 
    mutate(diff = -1*diff) %>% # invert so high-low
    filter(str_detect(coef, "jobs")) %>% 
    mutate(subgroup = "Containment policies"),
  
  gender_mods$model %>% 
    wald_test_submodel(.) %>% 
    mutate(diff = -1*diff) %>% # invert so Women - Male
    filter(str_detect(coef, "lockdown")) %>% 
    mutate(subgroup = "Gender"),
  
  ideology$model %>% 
    wald_test_submodel(.) %>% 
    filter(str_detect(coef, "lockdown")) %>% 
    mutate(subgroup = "Ideology"),
  
  vacc_refuse_mods$model %>% 
    wald_test_submodel(.) %>% 
    filter(str_detect(coef, "lockdown")) %>% 
    mutate(diff = -1*diff) %>% # invert so Refusers - Takers
    mutate(subgroup = "Vaccine refusal")
) %>% 
  
  mutate(att = case_when(str_starts(pattern = "gdp", string = coef) ~ "GDP growth",
                         str_starts(pattern = "jobs", string = coef) ~ "Job growth",
                         str_starts(pattern = "supplies", string = coef) ~ "Vaccine procurement",
                         str_starts(pattern = "lockdown", string = coef) ~ "Lockdown length",
                         str_starts(pattern = "deaths", string = coef) ~ "Deaths",
                         str_starts(pattern = "vaccinated", string = coef) ~ "Vaccination rate"),
         coef = case_when(str_starts(pattern = "gdp", string = coef) ~ gsub("gdp","GDP growth: ", coef),
                          str_starts(pattern = "jobs", string = coef) ~ gsub("jobs","Job growth: ", coef),
                          str_starts(pattern = "supplies", string = coef) ~ gsub("supplies","Vaccine procurement: ", coef),
                          str_starts(pattern = "lockdown", string = coef) ~ gsub("lockdown","Lockdown length: ", coef),
                          str_starts(pattern = "deaths", string = coef) ~ gsub("deaths","Deaths", coef),
                          str_starts(pattern = "vaccinated", string = coef) ~ gsub("vaccinated","Vaccination rate: ", coef))) %>% 
  mutate(coef = factor(coef, levels = rev(c("Job growth: +10%","Job growth: +5%","Job growth: -5%","Job growth: -10%",
                                            "Vaccination rate: 75%","Vaccination rate: 50%","Vaccination rate: 25%","Vaccination rate: 15%",
                                            "Lockdown length: 40 weeks","Lockdown length: 30 weeks","Lockdown length: 20 weeks")))) %>% 
  group_by(subgroup) %>% 
  mutate(lowest = diff == min(diff)) %>% 
  ungroup() %>% 
  mutate(subgroup = case_when(
    subgroup == "Containment policies" ~ "Containment policy: Stronger policies — Weaker policies",
    subgroup == "Food Poverty" ~ "Food poverty: More impacted — Less impacted ",
    subgroup == "Gender" ~ "Gender: Women — Men",
    subgroup == "Ideology" ~ "Ideology: Right-leaning — left-leaning",
    subgroup == "Vaccine refusal" ~ "Vaccine refusal: Refusers — Takers"
  )) 


ggplot(fig2_results, 
       aes(
         x = diff, 
         xmin = diff-1.96*diff_se, 
         xmax = diff+1.96*diff_se, 
         y = coef,
         color = subgroup,
       )) + 
  
  facet_wrap(subgroup ~ ., scales = "free_y", ncol = 1) +
  geom_point() +
  geom_errorbarh(height = 0.3) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_alpha_manual(values = c(0.3,1), guide = "none") +
  labs(title = "Difference in AMCEs conditioning on...", x = "Effect difference", y = "") +
  theme_bw() +
  theme(legend.position = "none",
        panel.background = element_blank())

ggsave("figures/figure2.pdf", width = 5.5, height = 7.5)

##### PCA SUMMARIES #####

###### TABLE S8 -- PCA PROPORTION OF VARIANCE SUMMARY #####

pca_mods <- read_rds("models/pca_recoding_mods.RDS")

max_components <- max(sapply(pca_mods, function (x) ncol(x$loadings)))

pca_mod_components <- data.frame(component = 1:max_components)

for (mod in names(pca_mods)) {
  pca_mod_components[[mod]] <- c(
    paste0(round(pca_mods[[mod]]$sdev^2/sum(pca_mods[[mod]]$sdev^2)*100,2),"%"),
    rep("", max_components - length(pca_mods[[mod]]$sdev))
  )
}

pca_mod_components %>% rename(
  "Food poverty" = "food",
  "Health score" = "who",
  "Containment policies" = "C",
  "Economic support policies" = "E",
  "Health compliance" = "health_compl",
  "Research spending" = "health_spend",
  "Component" = "component"
) %>% 
  xtable() %>% 
  print(only.contents = TRUE,
        include.rownames = FALSE,
        hline.after = c(0),
        booktabs = TRUE,
        file = "tables/pca_factors.tex")


###### TABLE S9-14 -- LOADINGS FOR EACH PCA MODEL #####
pca_mod_loadings <- lapply(pca_mods, loadings)

for (i in 1:length(pca_mod_loadings)) {
  
  pca_mod_loadings[[i]][,] %>% 
    as.data.frame() %>% 
    add_column(`Survey item` = rownames(.), .before = 1) %>%
    rename_with(function (x) gsub("Comp.","",x),starts_with("Comp")) %>% 
    xtable(digits = 3) %>% 
    print(only.contents = TRUE,
          booktabs = TRUE,
          hline.after = c(0), 
          include.rownames = FALSE,
          file = paste0("tables/pca_loadings_",names(pca_mod_loadings[i]),".tex"))
  
}

##### FIGURE S30 -- CHECK PLAUSIBILITY OF NO CARRYOVER EFFECT #####

round_mods <- lapply(1:8, function (r) {
  
  conjoint_lm(data = conjoint_data %>% filter(round == r), 
              country_fe = TRUE, round_fe = FALSE, 
              controls = mod_ctrls) %>% 
    format_model(., ref_levels, label = r)
  
}) %>% 
  do.call("rbind",.)

round_plot <- round_mods %>% 
  mutate(coef_full = paste0(att,": ",coef)) %>% 
  filter(non_ref == 1) %>% 
  ggplot(., aes(x = model, y = est, color = att, fill = att, group = coef_full,
                ymin = est-1.96*std_error, ymax = est+1.96*std_error)) +
  facet_grid(att ~ ., space = "free") +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_point(alpha = 0.8) +
  geom_line() +
  geom_ribbon(alpha = 0.15, color = NA) +
  labs(x = "Round", y = "AMCE estimate", color = "Attribute", fill = "Attribute") +
  theme_bw() +
  theme(legend.position = "bottom",
        text = element_text(size = 11),
        legend.text = element_text(size = 11),
        strip.text = element_text(size = 11))
plot(round_plot)
ggsave(round_plot, filename = "figures/round_plot.pdf", width = 7, height = 11)  


#### BENCHMARKING ANALYSIS ####

## NET OPINIONS

subj_pos <- conjoint_data %>% 
  mutate(gdp_pos = case_when(gdp == "+10%" ~ 0,
                             gdp == "+5%" ~ -0.25,
                             gdp == "0%" ~ -.5,
                             gdp == "-5%" ~ -0.75,
                             gdp == "-10%" ~ -1),
         
         jobs_pos = case_when(jobs =="+10%" ~ 0,
                              jobs =="+5%" ~ -0.25,
                              jobs == "0%" ~ -0.5,
                              jobs == "-5%" ~ -0.75,
                              jobs == "-10%" ~ -1),
         
         supplies_pos = ifelse(str_starts(supplies, "Quick"), 0, -1),
         
         lockdown_pos = case_when(lockdown =="10 weeks" ~ 0,
                                  lockdown =="20 weeks" ~ -0.333,
                                  lockdown =="30 weeks" ~ -0.666,
                                  lockdown =="40 weeks" ~ -1),
         
         vaccinated_pos = case_when(vaccinated == "75%" ~ 0,
                                    vaccinated == "50%" ~ -0.25,
                                    vaccinated == "25%" ~ -0.5,
                                    vaccinated == "15%" ~ -0.75,
                                    vaccinated == "5%" ~ -1),
         
         deaths_pos = case_when(deaths == "10 per million" ~ 0,
                                deaths == "30 per million" ~ -0.25,
                                deaths == "50 per million" ~ -0.5,
                                deaths == "70 per million" ~ -0.75,
                                deaths == "90 per million" ~ -1)) %>%
  
  mutate(profile_pos = gdp_pos + jobs_pos + supplies_pos + lockdown_pos + vaccinated_pos + deaths_pos, na.rm = TRUE) %>% 
  group_by(c_id, country) %>% 
  summarise(subj_pos = mean(profile_pos, na.rm = TRUE),
            gdp_pos = mean(gdp_pos, na.rm = TRUE),
            jobs_pos = mean(jobs_pos, na.rm = TRUE),
            supplies_pos = mean(supplies_pos, na.rm = TRUE),
            lockdown_pos = mean(lockdown_pos, na.rm = TRUE),
            vaccinated_pos = mean(vaccinated_pos, na.rm = TRUE),
            deaths_pos = mean(deaths_pos, na.rm = TRUE))

## Difference between home and overseas

home_overseas_diff <- conjoint_data %>% 
  left_join(subj_pos, by = c("country","c_id")) %>% 
  select(c_id, country, starts_with("politics_"), gov_relect, 
         all_of(mod_ctrls), 
         ideology, POLCONV_VDEM, system, ends_with("_pos"),
         health_spend_pca, health_compl_pca, intl_flights,
         POLCONV_VDEM,  system) %>% 
  mutate(marital_status = case_when(marital_status == "Yes" ~ 1,
                                    marital_status == "No" ~ 0,
                                    TRUE ~ NA),
         dep_children = case_when(dep_children == "Yes" ~ 1,
                                  dep_children == "No" ~ 0,
                                  TRUE ~ NA)) %>% 
  distinct()

for (i in seq(1,9,2)) {
  home_overseas_diff[paste0("home_bias_",as.integer((i+1)/2))] <- 
    home_overseas_diff[paste0("politics_",i)] - home_overseas_diff[paste0("politics_",i+1)]
}

home_overseas_long <- home_overseas_diff %>%  
  select(starts_with("home_bias")) %>% 
  pivot_longer(everything())

###### FIGURE 6 #####

country_net <- home_overseas_diff %>%  
  select(country, starts_with("home_bias")) %>% 
  pivot_longer(starts_with("home_bias")) %>% 
  group_by(name, country) %>%  
  summarise(diff = mean(value, na.rm = TRUE),
            se = sd(value, na.rm = TRUE)/sqrt(sum(!is.na(value))),
            non_missing = sum(!is.na(value)),
            n = n()) %>% 
  mutate(type = case_when(str_detect(name, "_1") ~ "Handling of Covid-19 pandemic",
                          str_detect(name, "_2") ~ "Handling of lockdown policy",
                          str_detect(name, "_3") ~ "Handling of Covid-19 vaccine campaign",
                          str_detect(name, "_4") ~ "Handling of economic policies",
                          str_detect(name, "_5") ~ "Handling of Covid-19 related deaths")) %>% 
  mutate(region = case_when(country %in% c("GHA", "UGA", "ZAF") ~ "Africa",
                            country %in% c("FR","IT","UK","SP") ~ "Europe",
                            country %in% c("AUS","CHN","IND", "JPN") ~ "APac",
                            country %in% c("BR","CAN","CHL","COL","US") ~ "Americas")) %>% 
  
  
  mutate(type_short = case_when(str_detect(name, "_1") ~ "Pandemic",
                                str_detect(name, "_2") ~ "Lockdown",
                                str_detect(name, "_3") ~ "Vaccines",
                                str_detect(name, "_4") ~ "Economy",
                                str_detect(name, "_5") ~ "Deaths"),
         name = "",
         country = factor(country,
                          levels = c("GHA","UGA","ZAF",
                                     "BR","CAN","CHL","COL","US",
                                     "AUS","CHN","IND","JPN",
                                     "FR","IT","SP","UK")))

base_plot <- ggplot(country_net %>% arrange(country), aes(y = type_short, x = diff, color = country)) +
  geom_point(position = position_dodge(-0.5)) +
  facet_grid(type_short ~ region, scales = "free") +
  geom_errorbarh(aes(xmin = diff-1.96*se, xmax = diff+1.96*se),
                 height = 0.2,
                 position = position_dodge(-0.5)) +
  geom_vline(xintercept = 0, linetype="dashed") +
  labs(x = "Net difference in opinion of own vs other governments",
       y = "",
       color = "") +
  theme_bw() +
  scale_color_manual(values = plot_colours) +
  theme(legend.position = "none",
        strip.text.y = element_blank())

regions <- c("Africa","Americas","APac","Europe")

region_plots <- list()
for (region_str in regions) {
  region_plots[[region_str]] <- country_net %>% 
    filter(region == region_str) %>%
    ggplot(aes(y = type_short, x = diff, color = country)) +
    geom_point() +
    scale_color_manual(values = plot_colours, name = "") +
    labs(color = region_str) +
    theme_bw() +
    theme(legend.text = element_text(size = 11),
          legend.background = element_rect(fill = "transparent"))
}

plot_grid(
  base_plot,
  plot_grid(
    get_legend(region_plots[["Africa"]]),
    get_legend(region_plots[["Americas"]]),
    get_legend(region_plots[["APac"]]),
    get_legend(region_plots[["Europe"]]),
    nrow = 1, ncol = 4,
    axis = "t",
    align = "hv",
    rel_widths = c(1.5,1,1,1)
  ), 
  nrow = 2, 
  rel_heights = c(5,2)
)

ggsave("figures/figure6.pdf", width = 8.25, height = 4.6)

###### FIGURE S32 -- BY IDEOLOGY ######

home_overseas_diff %>% 
  mutate(ideology_bin = ifelse(ideology <= 5, "Left/Centre","Right")) %>% 
  net_diff_subset(., "ideology_bin", "Ideology", plot_colours)

ggsave("figures/net_difference_by_ideology.pdf", width = 8, height = 7)

###### FIGURE S33 -- BY EDUCATION ######

home_overseas_diff %>% 
  filter(EDUCATION_LEVEL != "Unknown/missing") %>% 
  group_by(country, EDUCATION_LEVEL) %>% 
  nest() %>% 
  mutate(N = map(data, ~nrow(.x))) %>% 
  unnest(N) %>% 
  filter(N > 20) %>% 
  unnest(data) %>% 
  ungroup() %>% 
  mutate(EDUCATION_LEVEL = case_when(EDUCATION_LEVEL == "Less than primary completed" ~ "Less than primary",
                                     EDUCATION_LEVEL == "Primary completed" ~ "Primary",
                                     EDUCATION_LEVEL == "Secondary completed" ~ "Secondary",
                                     EDUCATION_LEVEL == "University completed" ~ "University")) %>% 
  net_diff_subset(., "EDUCATION_LEVEL", "Education level", plot_colours)

ggsave("figures/net_difference_by_education.pdf", width = 8, height = 7)

###### FIGURE S34 -- BY INCOME ######

home_overseas_diff %>% 
  mutate(income_abovemed = ifelse(income_abovemed == 1, "Yes", "No")) %>% 
  net_diff_subset(., "income_abovemed", "Income above median?", plot_colours)

ggsave("figures/net_difference_by_income.pdf", width = 8, height = 7)

###### FIGURE S35 -- BY HEALTH MEASURE COMPLIANCE ######

net_diff_subset <- function(data, subset_var, slab, plot_colours = plot_colours) {
  
  
  data[["subset"]] <- data[[subset_var]]
  
  calc_diff <- data %>%  
    select(country, starts_with("home_bias"),subset) %>% 
    filter(!is.na(subset)) %>% 
    pivot_longer(starts_with("home_bias")) %>% 
    group_by(name, country, subset) %>%  
    summarise(diff = mean(value, na.rm = TRUE),
              se = sd(value, na.rm = TRUE)/sqrt(sum(!is.na(value))),
              non_missing = sum(!is.na(value)),
              n = n()) %>% 
    mutate(type = case_when(str_detect(name, "_1") ~ "Handling of Covid-19 pandemic",
                            str_detect(name, "_2") ~ "Handling of lockdown policy",
                            str_detect(name, "_3") ~ "Handling of Covid-19 vaccine campaign",
                            str_detect(name, "_4") ~ "Handling of economic policies",
                            str_detect(name, "_5") ~ "Handling of Covid-19 related deaths")) %>% 
    mutate(region = case_when(country %in% c("GHA", "UGA", "ZAF") ~ "Africa",
                              country %in% c("FR","IT","UK","SP") ~ "Europe",
                              country %in% c("AUS","CHN","IND", "JPN") ~ "APac",
                              country %in% c("BR","CAN","CHL","COL","US") ~ "Americas")) %>% 
    
    
    mutate(type_short = case_when(str_detect(name, "_1") ~ "Pandemic",
                                  str_detect(name, "_2") ~ "Lockdown",
                                  str_detect(name, "_3") ~ "Vaccines",
                                  str_detect(name, "_4") ~ "Economy",
                                  str_detect(name, "_5") ~ "Deaths"),
           name = "",
           country = factor(country,
                            levels = c("GHA","UGA","ZAF",
                                       "BR","CAN","CHL","COL","US",
                                       "AUS","CHN","IND","JPN",
                                       "FR","IT","SP","UK")))
  
  base_plot <- ggplot(calc_diff %>% arrange(country), aes(y = type_short, x = diff, color = country, shape = subset)) +
    geom_point(position = position_dodge(-0.5)) +
    facet_grid( ~ region, scales = "free") +
    geom_errorbarh(aes(xmin = diff-1.96*se, xmax = diff+1.96*se),
                   height = 0.2,
                   position = position_dodge(-0.5)) +
    geom_vline(xintercept = 0, linetype="dashed") +
    labs(x = "Net difference in opinion of own vs other governments",
         y = "",
         color = "",
         shape = slab) +
    scale_color_manual(values = plot_colours, guide = "none") + 
    guides(shape = guide_legend(override.aes=list(color = "grey50"))) +
    theme_bw() +
    theme(legend.position = "top",
          legend.text = element_text(size = 11),
          legend.box.background = element_rect(color = "grey50", size = 2))
  
  regions <- c("Africa","Americas","APac","Europe")
  
  region_plots <- list()
  for (region_str in regions) {
    region_plots[[region_str]] <- calc_diff %>% 
      filter(region == region_str) %>%
      ggplot(aes(y = type_short, x = diff, color = country)) +
      geom_point() +
      scale_color_manual(values = plot_colours, name = "") +
      labs(color = region_str) + 
      theme(text = element_text(size = 11),
            legend.text = element_text(size = 11),
            strip.text = element_text(size = 11))
  }
  
  plot_grid(
    base_plot,
    plot_grid(
      get_legend(region_plots[["Africa"]]),
      get_legend(region_plots[["Americas"]]),
      get_legend(region_plots[["APac"]]),
      get_legend(region_plots[["Europe"]]),
      nrow = 1, ncol = 4,
      axis = "t",
      align = "hv",
      rel_widths = c(1.5,1,1,1)
    ), 
    nrow = 2, 
    rel_heights = c(7,2)
  )
  
}

home_overseas_diff %>% 
  mutate(health_compl_bin = ifelse(health_compl_pca <= median(health_compl_pca, na.rm = TRUE), "Low","High")) %>% 
  net_diff_subset(., "health_compl_bin", "Health measure compliance", plot_colours)

ggsave("figures/net_difference_by_health_compl.pdf", width = 8, height = 7)


###### FIGURE S36 -- BY HEALTH RESEARCH ATTITUDES ######

home_overseas_diff %>% 
  mutate(health_spend_bin = ifelse(health_spend_pca <= median(health_spend_pca, na.rm = TRUE), "Low","High")) %>% 
  net_diff_subset(., "health_spend_bin", "Health research spending attitudes", plot_colours)

ggsave("figures/net_difference_by_health_spend.pdf", width = 8, height = 7)

###### FIGURE S37 -- BY NO. INTERNATIONAL FLIGHTS ######

home_overseas_diff %>% 
  mutate(intl_flights = case_when(intl_flights == 0 ~ "None",
                                  intl_flights < 3 ~ "1-2",
                                  intl_flights < 5 ~ "3-4",
                                  intl_flights >= 5 ~ ">5")) %>% 
  mutate(intl_flights = factor(intl_flights,
                               levels = c("None","1-2","3-4",">5"))) %>% 
  group_by(country, intl_flights) %>% 
  nest() %>% 
  mutate(N = map(data, ~nrow(.x))) %>% 
  unnest(N) %>% 
  filter(N > 20) %>% 
  unnest(data) %>% 
  ungroup() %>% 
  net_diff_subset(., "intl_flights", "No. of international flights", plot_colours)

ggsave("figures/net_difference_by_intl_flights.pdf", width = 8, height = 7)

###### FIGURE S38 -- BY SUBJECT POSITIVITY ######

pos_net <- home_overseas_diff %>% 
  mutate(subj_pos_bin = ifelse(subj_pos < median(subj_pos), "Low", "High")) %>% 
  pivot_longer(starts_with("home_bias")) %>% 
  group_by(name, subj_pos_bin) %>%  
  summarise(diff = mean(value, na.rm = TRUE),
            se = sd(value, na.rm = TRUE)/sqrt(sum(!is.na(value))),
            non_missing = sum(!is.na(value)),
            n = n()) %>% 
  mutate(type = case_when(str_detect(name, "_1") ~ "Handling of Covid-19 pandemic",
                          str_detect(name, "_2") ~ "Handling of lockdown policy",
                          str_detect(name, "_3") ~ "Handling of Covid-19 vaccine campaign",
                          str_detect(name, "_4") ~ "Handling of economic policies",
                          str_detect(name, "_5") ~ "Handling of Covid-19 related deaths")) %>% 
  
  
  mutate(type_short = case_when(str_detect(name, "_1") ~ "Pandemic",
                                str_detect(name, "_2") ~ "Lockdown",
                                str_detect(name, "_3") ~ "Vaccines",
                                str_detect(name, "_4") ~ "Economy",
                                str_detect(name, "_5") ~ "Deaths"),
         name = "")

ggplot(pos_net %>% arrange(subj_pos_bin), aes(y = type_short, x = diff, color = subj_pos_bin)) +
  geom_point(position = position_dodge(-0.5), size = 2) +
  geom_errorbarh(aes(xmin = diff-1.96*se, xmax = diff+1.96*se),
                 height = 0.2,
                 position = position_dodge(-0.5)) +
  geom_vline(xintercept = 0, linetype="dashed") +
  labs(x = "Net difference in opinion of own vs other governments",
       y = "",
       color = "Positivity of conjoint treatments") +
  theme_bw() +
  scale_color_manual(values = c("firebrick2", "cadetblue4")) +
  plot_theme +
  theme(legend.position = "bottom")

ggsave("figures/net_difference_by_pos.pdf", width = 8, height = 6.5)

##### FIGURE S39 -- BY SUBJECT POSITIVITY, WITHIN COUNTRY #####

home_overseas_diff %>% 
  mutate(subj_pos_bin = ifelse(subj_pos < median(subj_pos), "Low", "High")) %>% 
  net_diff_subset(., "subj_pos_bin", "Positivity of conjoint treatments", plot_colours)

ggsave("figures/net_difference_by_pos_country.pdf", width = 8, height = 6.5)


###### FIGURE S40 -- BY REGIME TYPE ######

system_net <- home_overseas_diff %>% 
  filter(country != "FR") %>% 
  pivot_longer(starts_with("home_bias")) %>% 
  group_by(name, system) %>%  
  summarise(diff = mean(value, na.rm = TRUE),
            se = sd(value, na.rm = TRUE)/sqrt(sum(!is.na(value))),
            non_missing = sum(!is.na(value)),
            n = n()) %>% 
  mutate(type = case_when(str_detect(name, "_1") ~ "Handling of Covid-19 pandemic",
                          str_detect(name, "_2") ~ "Handling of lockdown policy",
                          str_detect(name, "_3") ~ "Handling of Covid-19 vaccine campaign",
                          str_detect(name, "_4") ~ "Handling of economic policies",
                          str_detect(name, "_5") ~ "Handling of Covid-19 related deaths")) %>% 
  
  
  mutate(type_short = case_when(str_detect(name, "_1") ~ "Pandemic",
                                str_detect(name, "_2") ~ "Lockdown",
                                str_detect(name, "_3") ~ "Vaccines",
                                str_detect(name, "_4") ~ "Economy",
                                str_detect(name, "_5") ~ "Deaths"),
         name = "")

ggplot(system_net %>% arrange(system), aes(y = type_short, x = diff, color = system)) +
  geom_point(position = position_dodge(-0.5), size = 2) +
  geom_errorbarh(aes(xmin = diff-1.96*se, xmax = diff+1.96*se),
                 height = 0.2,
                 position = position_dodge(-0.5)) +
  geom_vline(xintercept = 0, linetype="dashed") +
  labs(x = "Net difference in opinion of own vs other governments",
       y = "",
       color = "") +
  scale_color_manual(values = c("darkorange1", "darkorchid3", "darkolivegreen")) +
  plot_theme + 
  theme(legend.position = "bottom")


ggsave("figures/net_difference_by_regime.pdf", width = 8, height = 6.5)

###### TABLE S29 -- REGRESSION OF NET APPROVAL ON DEMOGRAPHICS #####


net_mods <- list()
net_fs <- list()
for (i in 1:5) {
  outc <- paste0("home_bias_",i)
  
  net_mod_data <- home_overseas_diff %>% 
    select(
      all_of(c(outc, "country",
               "age", "ideology", "gender", "dep_children", 
               "marital_status", "EDUCATION_LEVEL",
               "intl_flights", "health_compl_pca", 
               "health_spend_pca", "income_abovemed"))) %>% 
    drop_na()
  
  full_mod <- lm(paste0(outc, " ~ ."), 
                 data = net_mod_data)
  
  fe_mod <- lm(paste0(outc, " ~ country"), 
               data = net_mod_data)
  
  f_test <- anova(full_mod, fe_mod)
  
  net_mods[[i]] <- full_mod
  net_fs[[i]] <- c(f_test$F[2],f_test$`Pr(>F)`[2])
  
}

texreg::texreg(
  net_mods,
  table = FALSE, tabular = TRUE,
  custom.coef.map = list("age" = "Age",
                         "ideology" = "Ideology",
                         "genderMale" = "Gender: male",
                         "genderOther" = "Gender: other",
                         "Prefer not to say" = "Gender: prefer not to say",
                         "dep_children" = "Dependent children",
                         "marital_status" = "Martial status",
                         "EDUCATION_LEVELPrimary completed" = "Primary education",
                         "EDUCATION_LEVELSecondary completed" = "Secondary education",
                         "EDUCATION_LEVELUniversity completed" = "University education",
                         "EDUCATION_LEVELUnkown/missing completed" = "Unkown/missing education",
                         "income_abovemed" = "Income above median",
                         "intl_flights" = "No. of international flights",
                         "health_compl_pca" = "Health compliance (PCA)",
                         "health_spend_pca" = "Health spending (PCA)"),
  digits = 3,
  custom.header = list("\\textbf{Evaluation Dimension}" = 1:5),
  custom.model.names = c("Pandemic", "Lockdown", "Vaccines","Economy","Deaths"),
  custom.gof.rows = list("Country FE?" = rep("Yes", 5)),
  booktabs = TRUE,
  use.packages = FALSE,
  include.rsquared = FALSE,
  file = "tables/net_approval_on_demographics.tex")


##### TABLE 1 -- APPROVALS ON REELECT #####

covars_of_interest <- c("ideology", "health_compl_pca","health_spend_pca","intl_flights")

reelect_df <- home_overseas_diff %>% 
  mutate(reelect = ifelse(gov_relect == "Yes",1,0)) %>% 
  select(
    matches("politics_[13579]$"), 
    starts_with("home_bias_"), 
    all_of(c("reelect","country")),
    all_of(mod_ctrls),
    all_of(covars_of_interest)) # %>% 
  # drop_na(matches("politics_|home_bias_|other_"))

for (i in seq(1,9,2)) {
  reelect_df[[paste0("other_",i)]] <- reelect_df[[paste0("politics_",i)]] - reelect_df[[paste0("home_bias_",(i+1)/2)]]
}

incumb_dom <- lm(reelect ~ politics_1 + politics_3 + politics_5 + politics_7 + politics_9 +
                   gender + age + EDUCATION_LEVEL + marital_status + dep_children + income_abovemed + country,
                 reelect_df %>% mutate(across(starts_with("politics_"), function (x) x/10)))

incumb_net <- lm(reelect ~ home_bias_1 + home_bias_2 + home_bias_3 + home_bias_4 + home_bias_5 +
                   gender + age + EDUCATION_LEVEL + marital_status + dep_children + income_abovemed + country,
                 reelect_df %>% mutate(across(starts_with("home_bias_"), function (x) x/10)))

incumb_both <- lm(reelect ~ politics_1 + politics_3 + politics_5 + politics_7 + politics_9 +
                    other_1 + other_3 + other_5 + other_7 + other_9 +
                    gender + age + EDUCATION_LEVEL + marital_status + dep_children + income_abovemed + country,
                  reelect_df %>% mutate(across(starts_with(c("politics_", "other_")), function (x) x/10)))


texreg::texreg(list(incumb_dom, incumb_net, incumb_both), 
               table = FALSE, tabular = TRUE,
               custom.coef.map = list("politics_1" = "Handling of Covid-19 pandemic",
                                      "politics_3" = "Handling of lockdown policy",
                                      "politics_5" = "Handling of Covid-19 vaccine campaign",
                                      "politics_7" = "Handling of economic policies",
                                      "politics_9" = "Handling of Covid-19 related deaths",
                                      "home_bias_1" = "Handling of Covid-19 pandemic",
                                      "home_bias_2" = "Handling of lockdown policy",
                                      "home_bias_3" = "Handling of Covid-19 vaccine campaign",
                                      "home_bias_4" = "Handling of economic policies",
                                      "home_bias_5" = "Handling of Covid-19 related deaths",
                                      "other_1" = "Covid-19 pandemic (other governments)",
                                      "other_3" = "Lockdown policy (other governments)",
                                      "other_5" = "Covid-19 vaccine campaign (other governments)",
                                      "other_7" = "Economic policies (other governments)",
                                      "other_9" = "Covid-19 related deaths (other governments)"),
               digits = 3,
               custom.header = list("\\textbf{Evaluation}" = 1:3),
               custom.model.names = c("\\textit{Domestic}","\\textit{Net}","\\textit{Domestic + Other}"),
               custom.gof.rows = list("Subject controls?" = rep("Yes", 3),
                                      "Country FE?" = rep("Yes", 3)),
               booktabs = TRUE,
               use.packages = FALSE,
               include.rsquared = FALSE,
               file = "tables/incumb_approval_regs.tex")


###### TABLE S21 -- WALD TEST #####

wald_comps <- data.frame(
  var = names(incumb_dom$coefficients)[2:6],
  dom_net = wald_test(incumb_dom, incumb_net),
  dom_both = wald_test(incumb_dom, incumb_both),
  net_both = wald_test(incumb_net, incumb_both)
) %>% 
  mutate(var = case_when(var %in% c("politics_1","home_bias_1") ~ "Covid-19 pandemic",
                         var %in% c("politics_3","home_bias_2") ~ "Lockdown policy",
                         var %in% c("politics_5","home_bias_3") ~ "Vaccine campaign",
                         var %in% c("politics_7","home_bias_4") ~ "Economic policies",
                         var %in% c("politics_9","home_bias_5") ~ "Covid-19 related deaths")) %>% 
  mutate(across(where(is.numeric), function (x) round(x,3)),
         across(where(is.numeric), function (x) case_when(abs(x) > 3.291 ~ paste0(x,"***"),
                                                          abs(x) > 2.576 ~ paste0(x,"**"),
                                                          abs(x) > 1.960 ~ paste0(x,"*"),
                                                          TRUE ~ paste0(x)))) %>% 
  xtable(.) %>% 
  print(only.contents = TRUE,
        booktabs = TRUE,
        include.rownames = FALSE,
        include.colnames = FALSE,
        hline.after = NULL,
        file = "tables/wald_coefs.tex")

###### TABLE S22 -- VIF SCORE CHECK #####

vif_scores <- lapply(list(incumb_dom, incumb_net, incumb_both), function (x) {
  
  car::vif(x) %>% as.data.frame() %>% 
    mutate(variable = rownames(.), 
           mod = case_when("other_1" %in% variable ~ "Domestic + Other",
                           "home_bias_1" %in% variable ~ "Net",
                           "politics_1" %in% variable ~ "Domestic"),
           variable = case_when(variable %in% c("politics_1","home_bias_1") ~ "Covid-19 pandemic",
                                variable %in% c("politics_3","home_bias_2") ~ "Lockdown policy",
                                variable %in% c("politics_5","home_bias_3") ~ "Vaccine campaign",
                                variable %in% c("politics_7","home_bias_4") ~ "Economic policies",
                                variable %in% c("politics_9","home_bias_5") ~ "Covid-19 related deaths",
                                
                                variable %in% c("other_1") ~ "Covid-19 pandemic (other)",
                                variable %in% c("other_3") ~ "Lockdown policy (other)",
                                variable %in% c("other_5") ~ "Vaccine campaign (other)",
                                variable %in% c("other_7") ~ "Economic policies (other)",
                                variable %in% c("other_9") ~ "Covid-19 related deaths (other)",
                                TRUE ~ variable))
}) %>% 
  do.call("rbind",.) %>% 
  select(all_of(c("GVIF","variable","mod"))) %>% 
  pivot_wider(names_from = mod, values_from = GVIF) %>% 
  mutate(variable = str_to_sentence(variable),
         variable = case_when(variable == "Education_level" ~ "Education level",
                              variable == "Marital_status"~ "Marital status",
                              variable == "Dep_children" ~ "Dep. children",
                              variable == "Income_abovemed" ~ "Income above median",
                              TRUE ~ variable)) %>% 
  xtable(.) %>% 
  print(digits = 3, 
        only.contents = TRUE,
        include.rownames = FALSE,
        include.colnames = FALSE,
        hline.after = NULL,
        file = "tables/vif_scores.tex")

##### APPROVALS ON REELECT SUBSETS ####
###### TABLE S23 -- BY EDUCATION LEVEL #####

reelect_df %>% 
  filter(EDUCATION_LEVEL != "Unknown/missing") %>% 
  approval_subset(.,
                  "EDUCATION_LEVEL", 
                  "Reported compliance with health measures",
                  c("Univ.","Second.","Prim.","\\textless Prim."),
                  split_at_median = FALSE
  )

###### TABLE S24 -- BY IDEOLOGY ######

reelect_df %>% 
  approval_subset(.,
                  "ideology", 
                  "Ideology",
                  c("Right","Left/Centre"),
                  split_at_median = TRUE
  )

###### TABLE S25 -- BY INCOME ABOVE MEDIAN #####

reelect_df %>% 
  filter(!is.na(income_abovemed)) %>% 
  approval_subset(.,
                  "income_abovemed", 
                  "Income > median?",
                  c("Yes","No"),
                  split_at_median = FALSE
  )

###### TABLE S26 -- BY HEALTH COMPLIANCE #####

reelect_df %>% 
  approval_subset(.,
                  "health_compl_pca", 
                  "Reported compliance with health measures",
                  c("Low","High"),
                  split_at_median = TRUE
  )

###### TABLE S27 -- BY HEALTH SPENDING INTEREST #####

reelect_df %>% 
  approval_subset(.,
                  "health_spend_pca", 
                  "Agreement with greater health research spending",
                  c("High","Low"),
                  split_at_median = TRUE
  )

###### TABLE S28 -- BY NO. INTERNATIONAL FLIGHTS ######

reelect_df %>% 
  mutate(intl_flights = case_when(intl_flights == 0 ~ "None",
                                  intl_flights < 3 ~ "1-2",
                                  intl_flights < 5 ~ "3-4",
                                  intl_flights >= 5 ~ "\\textgreater 5")) %>% 
  mutate(intl_flights = factor(intl_flights,
                               levels = c("None","1-2","3-4","\\textgreater 5"))) %>% 
  arrange(intl_flights) %>% 
  approval_subset(.,
                  "intl_flights", 
                  "Number of international flights (2019)",
                  split_at_median = FALSE
  )



#### SUPPLEMENTARY ANALYSES ####
##### TABLE S1 -- ELECTION DATES #####
## Produced by hand 
##### TABLE S18 -- Shift in predicted probabilities #####
avg_net_shift <- reelect_df %>% 
  filter(country!="CHN") %>% 
  mutate(across(starts_with("home_bias_"), function (x) x/10)) %>% 
  select(all_of(c("country")), matches("home_bias_[1-5]")) %>% 
  pivot_longer(-country) %>% 
  group_by(name, country) %>% 
  summarise(avg = mean(value, na.rm = TRUE),
            sd = sd(value, na.rm = TRUE)) %>% 
  
  mutate(shift = incumb_net$coefficients[name]*sd) %>% 
  
  pivot_wider(id_cols = "name",
              names_from = "country",
              values_from = "shift") %>% 
  mutate(name = case_when(
    name == "home_bias_1" ~ "Covid-19 pandemic",
    name == "home_bias_2" ~ "Lockdown policy",
    name == "home_bias_3" ~ "Covid-19 vaccine campaign",
    name == "home_bias_4" ~ "Economic policies",
    name == "home_bias_5" ~ "Covid-19 related deaths")
  ) %>% 
  rename(" " = name)

xtable(avg_net_shift, digits = 3) %>% 
  print(include.rownames = FALSE,
        only.contents = TRUE,
        hline.after = 0,
        booktabs = TRUE,
        file = "tables/reelect_pred_probs.tex")

##### TABLE S19 -- Net approval bivariate models #####

reelect_net_lms <- lapply(1:5, function (i) {
  lm(paste0("reelect ~ home_bias_",i," + gender + age + EDUCATION_LEVEL + income_abovemed + dep_children + marital_status + country"), 
     reelect_df %>% mutate(across(starts_with("home_bias_"), function (x) x/10)))
}
)

reelect_net_lms[[6]] <- incumb_net

texreg::texreg(reelect_net_lms, table = FALSE, tabular = TRUE,
               custom.coef.map = list("home_bias_1" = "Handling of Covid-19 pandemic",
                                      "home_bias_2" = "Handling of lockdown policy",
                                      "home_bias_3" = "Handling of Covid-19 vaccine campaign",
                                      "home_bias_4" = "Handling of economic policies",
                                      "home_bias_5" = "Handling of Covid-19 related deaths",
                                      "Baseline re-election rate" = "Intercept"),
               digits = 3,
               custom.model.names = c("(1)","(2)","(3)","(4)","(5)","(6)"),
               custom.gof.rows = list("Subject controls?" = rep("Yes", 6),
                                      "Country FE?" = rep("Yes", 6)),
               booktabs = TRUE,
               use.packages = FALSE,
               file = "tables/reelection_regressions_net.tex")

##### TABLE S20 -- Net results by country #####
reelect_country_net <- reelect_df %>% 
  mutate(across(starts_with("home_bias_"), function (x) x/10)) %>% 
  group_by(country) %>% 
  filter(country != "CHN") %>% 
  nest() %>% 
  mutate(mod = map(data, ~lm(reelect ~ home_bias_1 + home_bias_2 + home_bias_3 + home_bias_4 + home_bias_5 +
                               gender + age + EDUCATION_LEVEL + income_abovemed + dep_children + marital_status,
                             data = .)))

reelect_country_net_mods <- reelect_country_net %>% pull(mod)
names(reelect_country_net_mods) <- reelect_country_net$country

texreg::texreg(reelect_country_net_mods,
               sideways = TRUE, table = FALSE, tabular = TRUE,
               custom.coef.map = list("home_bias_1" = "Handling of Covid-19 pandemic",
                                      "home_bias_2" = "Handling of lockdown policy",
                                      "home_bias_3" = "Handling of Covid-19 vaccine campaign",
                                      "home_bias_4" = "Handling of economic policies",
                                      "home_bias_5" = "Handling of Covid-19 related deaths",
                                      "Baseline re-election rate" = "Intercept"),
               digits = 3,
               custom.gof.rows = list("Subject controls?" = rep("Yes", 15)),
               booktabs = TRUE,
               use.packages = FALSE,
               no.margin = TRUE,
               file = "tables/reelection_country_regressions_net.tex")


